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ABSTRACT 

Numerical simulations of galaxy formation require a number of parameters. Some of 
these are intrinsic to the numerical integration scheme (e.g. the timestep), while oth- 
ers describe the physical model (e.g. the gas metallicity) . In this paper, we present 
results of a systematic exploration of the effects of varying a subset of these param- 
eters on simulations of galaxy formation. We use iV-body and "Smoothed Particle 
Hydrodynamics" techniques to follow the evolution of cold dark matter and gas in a 
small volume. We compare a fiducial model to 24 different simulations, in which one 
parameter at a time is varied, focussing on properties such as the relative fraction of 
hot and cold gas, and the abundance and masses of galaxies. We find that for reason- 
able choices of numerical values, many parameters have relatively little effect on the 
galaxies, with the notable exception of the parameters that control the resolution of 
the simulation and the efficiency with which gas cools. 
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1 INTRODUCTION 



Over tliti pdhL LwL) decddeb iiumeiical simulaLiuiib licivti be- 
come the most commonly used technique for investigating 
the formation of structure in the universe. Originally, A''- 
body simulations were used to calcula te the clustering evo- 



fc White 1995: lAnninos fc Norman 1996| : [Bryan fc Norman 
199^; Eke, Navarro & Frenk 1998). Since the cooling time 



lution of collisionless dark matter 



1985 



Centrella fc Melott 1985; Davis et al. 1985). Subse 
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quently, countless A-body studies have addressed problems 
ranging from the large-scale distribution of dark matter to 
the inner structure of dark matter haloes in a programme 
whose results underpin much of modern cosmology. 

The success of the A-body approach stimulated the ex- 
tension of simulation techniques to processes involving gas. 
Eulerian, Lagrangian and deformable grid methods have 
been developed for this purpose. One of the simplest cos- 
mological problem involving gas is th e evolution of the intr- 
acluster medium in rich clusters (e.g.lEvrard 1990|; [f homad 
fc Couchman 1992| |Cen fc Ostriker 1994^ [Navarro, Frenk 



of the gas is longer than the Hubble time in all but the cen- 
tral regions of clusters, the intracluster gas may be approxi- 
mated as non-radiative. This is an important simplification 
that allows reliable and repeatable results to be obtained 
with a variety of techniques and numerical resolutions, as 
demonstrated by a recent systematic comparison based on a 
simulation o f a galaxy cluster in the cold dark matter (CDM) 
cosmogony (Frenk et al. 1999). 



Another cosmological problem that can be addressed 
using gas dynamics is the evolution of cold gas clouds at 
high-redshift. These are responsible for the Lyman-a forest 
lines seen in quasar spectra and their low temperatures and 
moderate density contrasts ensure that cooling processes are 
relatively unimportant. Simulations of the Lyman-a forest 
in CDM models have been very successful and have radically 
changed the prevaili ng view of t he high redshift gaseous uni- 
verse (e.g. Cen et al. 



1994 



1995; Hernquist et al. 1996; Katz 



Hernquist fc Weinberg 1999). Recent work, however, has 
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highlighted the need to take numerical effects carefully into 
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accoun t when making comparisons with observations (The- 
uns et ^l. 199^ ) 



Simulating the formation of galaxies requires the incor- 
poration of radiative coohng processes into the fluid equa- 
tions, enabling the gas to achieve high eno ugh density con- 
trasts so that stars can plausibly form (e .g. Katz, Hernquist 
& Weill ' 



berg 1992 



Cen & Ostrikcr 1992 



& Davi 5 1994 ; Pearce et al. 199S; ) . In general, the inclusion 



Evrard, Summers 



of such processes presents a formidable numerical challenge 
because in hierarchical clustering theories, such as CDM, 
cooling is most efficient in small objects at high redshift 
where cooling times are very short. As a result, the out- 
come of a simulation depends on its resolution. Fortunately, 
the problem is less intransigent than it might appear. It is 
clear that only a small fraction of gas can have cooled into 
small pre-galactic objects at high redshift. This is normally 
attributed to the effects of "feedback", a term frequently 
used to denote a generic mechanism that prevents t he gas 
from c ooling c atastrophicallv and turning into s tars (White 
& Rees| 1978|; |Cole 199l|; |White fc Frenk 199l|). Afihough 



feedback is poorly understood, it almost certainly involves 
stellar energy input (supernova explosions, galactic winds, 
etc. ) to the intergalactic medium, and, perhaps, even energy 
input from active galactic nuclei. 

At present, feedback effects need to be included in sim- 
ulations using a phenomenological model. In some cases, 
a prescription to turn some of the gas into stars and to 



199J 



S ;cinmctz 



Navarro & White 


cases (e.g. 


Pearce 



et al. 1!'99), the effect of resolution itself is used to cut off 



cooling; only gas in objects above the resolution limit of the 
simulation can cool efficiently. 

Although resolution is the most important numerical 
concern in simulations of galaxy formation, several other 
numerical and physical parameters influence the outcome of 
a simulation. Examples of the former are the gravitational 
softening or the size of the timestep; examples of the latter 
are the metallicity of the gas or the baryon fraction. Often 
parameter values are set in the absence of rigorous criteria. 
In this paper we undertake a systematic study of a limited 
subset of the parameters required in simulations of galaxy 
formation. We consider the way in which such parameters af- 
fect, for example, the amount of gas that cools, the number 
and masses of the galaxies that form, etc. This is not in- 
tended as an exhaustive investigation, but rather as a guide 
to the sensitivity of simulation results to these particular 
parameters. 

We restrict our study to one particular implementation 
of numerical hydrodynamics, t he "Smoothed Particle Hy- 
drodynarn ics" (SPH) technique ( Gingold fc Monaghan 1977 



Lucy 1977). This Lagrangian approach is well suited to stud- 
ies of galaxy formation because large dynamic variations in 
gas density and temperature are easily accommodated. In 
the implementation used in this paper, our code models pro- 
cesses such Eis adiabatic compression and expansion, shock 
heating and radiative cooling. For simplicity, we have cho- 
sen to ignore the effects of photo-ionization (studied pre- 
viously, for example, by Weinberg, Hernquist & Katz 1997 



and Navarro fc Steinmetz 1997j ) and feedback. We intend to 
include these effects in a subsequent set of tests. 

The rest of this paper is organized as follows. In Sec- 



tion 2 we outline the numerical methods used, particularly 
our implementation of gravity and hydrodynamics with ra- 
diative cooling. In Section 3 we define a "fiducial" simulation 
and describe its evolution, focussing on the distribution of 
gas in the temperature-density plane and on the properties 
of "galaxies." In Section 4 we compare this fiducial simula- 
tion to a set of simulations in which one parameter at a time 
is varied. In particular, we consider the correspondence of 
galaxies in different simulations as well as properties such 
as their mass distribution and the distribution of baryons 
within dark matter haloes. We summarise our results in Sec- 
tion 5 and conclude in Section 6. 



2 NUMERICAL METHOD 

The data analysed in this paper were generated using hydra 



P] ( |Couchman, Thomas fc Pearce 1995 , hereafter CTP95), an 
adaptive particle-particle/particle-mesh (ap'^m) code incor- 
porating SPH. The simulations were performed on single- 
processor workstations, although the same code has been 



implement ed in parallel on a Cray T3D (Pearce & Couch 



man 1997). Below we briefly outline how hydra works. 



2.1 Gravity 

The AP''m method for calculating the gravitational forces 
involves three distinct parts; the particle-mesh (PM), 
particle-particle (PP) and reflnement placing algorithms. 
The implementation used in hydra is detailed in Couch- 
man (1991), although we outline the key concepts below. 



The PM algorithm is an efficient method for calculat- 
ing forces accurately down to scales of the order of the 
Nyquist wavelength of the mesh used. The computational 
speed scales ~ 0(L'' logL), where is the total number of 
mesh cells. The density field is discretized by smoothing the 
particles on to a regular cubi c mesh using the Triangular - 
Shaped Cloud (TSC) kernel ( |Hockney fc Eastwood 198l|) . 



The potential of this distribution is then obtained by per- 
forming a Fast-Fourier Transform (FFT) , multiplying by the 
appropriate Green's function, and doing an inverse FFT. 
The potential at each mesh point is differenced to obtain 
the forces, which are then interpolated back to the particle 
positions using the same TSC kernel. 

The PM method is augmented on small scales by the 
PP method, which calculates pairwise forces directly via a 
neighbour search. This is performed using a coarse mesh 
and searching only the 27 nearest cells for neighbours. The 
coarser this mesh is relative to the FFT mesh, the more 
accurate the force calculation becomes at the expense of 
extra PP work. 

On very small scales, the gravitational force must be 
softened not only to avoid effects due to two-body interac- 
tions, but also to limit the need for very small timesteps. 
HYDRA uses a spline (S2) softening function which returns 
the Newtonian value of the force for r > s, where s is 
the softening length. For convenience we express this value 

t This code is in the public domain and can be obtained from 
the Hydra Consortium at the following URL : 
http: / /phobos. astro. uwo.ca/hydra_consort/ 
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as the more common Plummer softening, e (e.g. Binney 
Tremai le 1987), using the equivalence s = 2.34e. (This pro- 



duces asymptotically matching force laws at both large and 
small scales.) At late times we employ a fixed physical soft- 
ening which is the value normally quoted as the softening 
for the simulation. In comoving coordinates this will grow 
backwards in time proportional to {1 + z). To maintain force 
accuracy we do not allow the softening to grow indefinitely, 
introducing a maximum value, Smax, which is typically set 
to equal 0.3 FFT cells. Above a certain redshift (typically 
around z — 1) the softening is constant in comoving coordi- 
nates. 

Since gravity is an attractive force, the particle distri- 
bution evolves into a highly clustered state where large over- 
densities are inevitable. Hence the ratio of PP to PM work 
rapidly increases, leading to a substantial increase in the 
computatio nal effort for eac h step. The refinement placing 
algorithm (Couchman 1991) eases this by locating highly 
clustered regions and recursively placing submeshes (or re- 
finements) over these volumes. This technique drastically 
reduces the total amount of PP work per timestep, passing 
much of the effort over to the faster PM algorithm. 



2.2 Hydrodynamics 

Hydrodynamical forces are modelled using the SPH formal- 
ism, which because of its Lagrangian nature, is well suited to 
problems with a large dynamic range in density - a natural 
occurrence in cosmological simulations. Solving the hydro- 
dynamic equations involves estimating various fields such as 
the density and pressure gradient, from the discrete parti- 
cle distribution. SPH tackles this by smoothing the particles 
over a finite volume using a spherically-symmetric smooth- 
ing kernel. To cope with large density fiuctuations, the range 
of the kernel is set for each particle to be the radius of a 
sphe re containin g approximately A^sph neighbouring parti- 
cles (Wood 1981). The neighbour search is performed at the 
same time as the PP calculation and hence introduces no 
additional computational overhead. This procedure intro- 
duces a maximum SPH search length equivalent to the size 
of a PP chaining cell which leads to a minimum resolved 
density. For normal simulation parameters this limit is ap- 
proximately the mean cosmological density. 

Stable integration of the fiuid equations requires the in- 



troduction of an artificial viscosity term ( Monaghan & Gin- 
gold 198J). The implementation used in this study is based 
on the divergence of the fiow, rather than the more usual 
relative pairwise velocity of particles. 

Some excellent reviews of the SPH techn i que are avai l- 
able in the literature (e.g. Monaghan 1992; [Benz 1990 ); 



the implementation in hydra is summarised in Thomas 
& Couchman (1992, hereafter TC92; see also CTP95). The 



reader is also referred to recent work by Thacker et al. (199J 



who have investigated different implementations of the SPH 
force calculation, specifically the force symmetrization and 
artificial viscosity assumptions. Part of that study (Section 
3.5) utilised the same initial conditions we use here and ex- 
amined the eff'ects of different SPH implementations. Our 
im pleme nt at ion is equivalent to version 1 of Thacker et 
al. ( |199^ . 




og (T/K) 



Figure 1. The fits used in HYDRA for the normalised cooling func- 
tion, A(T). The three curves illustrate this for three values of the 
metallicity: Z/Zq = 0.0, 0.5, 1.0. The function is only plotted for 
temperatures down to 10* K since this is the limit where HYDRA 
switches ofi^ the cooling. (Physically, negligible thermal energy is 
lost below this value.) Overlaid on the Z = 0.5Zq line is the in- 
terpolated curve from the tabulated values given by Sutherland 
& Dopita (1993). 



2.3 Radiative Cooling 

Radiative cooling requires adding a sink to the energy equa- 
tion derived from the emissivity function. 



e{n,T)^n^AiT), 



(2.1) 



where n is the number density of atoms/ions and electrons, 
and A(r) is the (temperature dependent) normalised cool- 
ing function. Our version of hydra uses a cooling function 
made up from a series of power law fits to the opticall y thin 
radiative cooling code of Raymond, Cox & Smith (197f) 
as detailed by TC92. This fit contains contributions from 
bremsstrahlung and from line cooling due to hydrogen, he- 
lium and heavier elements with an assumed abundance Z 
times the solar value. The fit is relatively good above 10''K 
- below this temperature the cooling function drops precip- 
itously and we take it to be zero. Figure |^ illustrates these 
curves for metallicities Z/Zq = 0.0, 0.5, 1. 0. Th e Z = 0.5Zq 
cooling function of Sutherland & Dopita (1993) is overlaid. 



2.4 Numerical Integration 

The numerical integration algorithm is described in CTP95, 
who evaluated several algorithms and settled for a sim- 
ple PEC (Predict, Evaluate, Correct) scheme. This has 
the advantage of keeping storage to a minimum while also 
allowing arbitrary changes in the timestep, if the force 
changes abruptly in time (as it can do in the vicinity of 
shocks). HYDRA evaluates the timestep. At, such that At — 
K min(0.4At„, 0. 25 Ata, 0.0625 Atn), where At„ = e/umax, 
Ata = (e/fflmax)^''^ aud Atn = 1/H, with u,nax being the 
current largest speed, amax the current largest acceleration 
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and H the Hubble parameter. The value of k is a normalisa- 
tion constant typically set to unity for adiabatic simulations 
(TC92). A Courant-style condition is implicit in Atv, since 
the minimum SPH smoothing length for the gas is set by 
the softening. Note that a single timestep is used to advance 
all particles. In principle, a further gain in efficiency could 
be obtained by implementing individual timesteps for each 
particle. 

Although the cooling timescale is often shorter than the 
dynamical (or sound-crossing) timescale, we do not shrink 
the timestep to follow cooling explicitly. Instead, we apply 
the cooling condition, assuming a constant density, at the 
end of the timestep. Thermal energy is removed such that 



Gi—Aei 



■At, 



Pi 



(2.2) 



with subscript i denoting the i"* particle and e, its specific 
thermal energy. This results in stable and physically reason- 
able behaviour, even for relatively long timesteps. 



3 THE FIDUCIAL SIMULATION 
3.1 Initial Conditions 

Throughout this paper we examine the effects of varying var- 
ious numerical and physical parameters relative to a fiducial 
simulation. Our fiducial initial conditions were chosen to cor- 
respond to the standard cold dark matter model (SCDM). 
Thus we choose the cosmic density parameter, SI = 1, the 
cosmological constant, A = and the Hubble constant, 
h = n Althniiffh this model is now out of favour (fnr 



exampl e, its power spectrum cannot simultaneouslv fit both 
the large-scale COBE normalisation and the amplitude of 
small-scale galaxy clustering), it is a well studied model with 
moderately large power on all scales relevant to galaxy for- 
mation. In any case, the precise form of the initial conditions 
is not important for the tests carried out here. 

The initial mass distribution was generated as described 
by Efstathiou et al. (1985), by perturbing a uniform distri- 



bution of particles using the Zel'dovich approximation. The 
amplitudes of the Gaussian random field were set using Bond 
& Efstathiou's (1984) approximation to the CDM power 
spectrum. The normalisation was set so that erg, the rms 
of linear mass fluctuations in top-hat spheres of comoving 
radius, 8/i~^Mpc, matches the value derived from t he local 
abundance of rich cluste r s. We adopted erg = . 6 (White, 
Efstathiou & Frcnk 19931 IVianna & Liddle 19961 lEke, Cok 



& Frenk 1996). The size of the simulation box (in comoving 
coordinates) was fixed at 10/i~^Mpc, with a mesh resolu- 
tion of L = 64 (implying a cell-length of 156.25 /i~^kpc). 
Although the box is not big enough to follow the evolution 
of large scale structure accurately (waves on the scale of the 
box are not evolving linearly at the final time), this is not 
important since we are not attempting to compare our re- 
sults directly to observations - the volume is far too small 
for such purposes. 



t J?o = lOO/ikms-lMpc-l 



Table 1. The values of the parameters chosen for the fiducial 
simulation that were subsequently varied. N is the total number 
of particles, A^sph is the number of gas particles over which the 
SPH smoothing is carried out, L is the size of the FFT mesh, 
•5max is the maximum comoving size of the softening in units of 
FFT cells, EQ is the efl^ective Plummer softening, re is the timestep 
normalisation, Hi, is the baryon fraction, Z is the gas metallicity, 
Zi is the initial redshift and Ti is the initial gas temperature. The 
choice of these parameters is discussed in more detail in the text. 



Parameter 


Value 


N 


2 X 32^ 




32 


L 


64 




0.3 


eo 


10/i-lkpc 


K 


1 




0.06 


z 


O.5Z0 




24 


Ti 


lOOK 



Table |l| summarises the fiducial values of all the parame- 
ters that we study in this paper. We represent the mass using 
32"^ gas and 32"^ dark matter particles. The initial softening 
was set to 0.3 FFT cells (equivalent to a comoving Plummer 
softening of ~ 20 /i~^kpc), but after z ~ 1 it switched to a 
fixed physical softening of 10 /i~^kpc. The number of neigh- 
bours over which the SPH algorithm smooths was set to 
TVspH = 32. The simulation was run from an initial epoch, 
Zi — 24, to the present {z — 0), with the data output at 
regular intervals. 

Regarding the physical parameters, we chose to set the 
baryon fraction equal to Qh — 0.06, con sistent with the re- 
sults from nu cleosynthesis calculations (Copi, Schramm & 



Turner 1995). The global gas metallicity was set to Z = 
0.5Zq, similar to the value o btained from spectral obser- 
vations of intracluster gas (e.g. Mushotzky et al. 1996| ). The 
gas was given a cold start, with the initial temperature set 
to lOOK . 



3.2 The Baryon Phase Diagram at z = 

The two objects that we study in this paper are the dark 
matter haloes and the cold, dense clumps of gas that we de- 
fine as galaxies. All our comparisons are carried out at z = 0. 
The selection criteria used to catalogue the haloes and the 
galaxies are described in detail in section 3.4, but as motiva- 



tion for the galaxy selection procedure, we initially develop 
an understanding of the correspondence between the spa- 
tial distribution and the thermal evolution of the baryons. 
Figure ^ illustrates the temperature-density distribution of 
the gas ai z = 0, spanning 7 orders of magnitude in density 
and 9 orders of magnitude in temperature. We have split 
the distribution into three phases, which we label as the 
uncollapsed phase, the shocked phase and the galaxy phase. 
Accompanying this figure are three dot-plots, showing the 
corresponding spatial distribution of each phase. 

3.2.1 The uncollapsed phase 

We define the uncollapsed phase as the regime in which gas 
particles have T < IK and p < 10"^ (p). The boundaries 
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Figure 2. (top-left) The temperature -density phase diagram of all baryons in the fiducial siirmlation at the final output time (2 = 0). 
The density is displayed in units of the mean gas density and the temperature in Kelvin. Positions of the gas particles obeying the various 
selection criteria defined in the text are shown in the other three plots, specifically the uncollapsed, shocked and galaxy phases. It is 
evident that each phase shows a considerable difference in structure: the uncollapsed phase particles are homogeneously distributed, the 
shocked phase particles show large fluctuations in density and the galaxy phase is exclusively made up of tight clumps of particles. 



for this phase are chosen such that they enclose the parti- 
cles that he on the power law relation, which is an adiabat. 
These particles are still largely in free expansion and are, on 
average, cooling adiabatically with T oc (1 + z)^ . Spatially, 
the particles are homogeneously distributed, as a result of 
occupying a small range in density. The hard lower limit 
on this quantity (visible for all temperatures) is artificial, 
imposed by the algorithm due to the maximum distance 
it can search for A^sph neighbours. For the fiducial simula- 
tion, Pmjjj ~ 1.9 (p), clearly an undesirably large value. As 
a result, the low density gas is forced to unphysically high 
values (the gas should have p < (p)). The latest version of 
HYDRA now properly accounts for low density gas and we 
have checked that this makes no significant difference to the 
properties of the high density material. 



3.2.2 The shocked phase 

The shocked gas phase contains particles with p < 10'' (p) 
and T > IK as well as all particles with T > 10® K. Gas 
entered this phase due to the collapse and violent relaxation 
of gravitationally unstable regions, in which the bulk kinetic 
energy of the gas was transferred, via shocking, into heat. As 
a result, the spatial distribution of this phase exhibits large 
fluctuations in density: halo atmospheres lie in this regime 
and the shocked gas is largely absent from the voids. 



3.2.3 The galaxy phase 

The particles in the galcixy phase lie in the regime, p > 
10^ (p), and T < lO^K . The strong horizontal line at lO^K 
is due to the fact that the cooling rate drops to zero at 
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Figure 3. The evolution of the uncollapsed, shocked and galaxy 
phases, plotted as the fraction of gas particles in each phase (with 
respect to the total mass of gas particles), as a function of the 
expansion factor, a = (1 + z)^^. The shaded region marks the 
time before the start of the simulation (a = 0.04). Redshift is 
indicated along the top of the figure. 



this temperature. Particles above a temperature of 10^ K 
must lose a significant fraction of their thermal energy before 
reaching the main locus, a transition made possible by the 
large contributions from hydrogen and helium line cooling 
at these temperatures. The plot of their positions highlights 
the fact that cold, dense gas forms very tight clumps, with 
spatial extents of the order of the gravitational softening. 
To select "galaxies" , we look at the particles that lie in this 
region of the phase diagram. 



3.3 Evolution of the Baryon Phase Diagram 

To illustrate the global evolution that led to the baryon 
phase diagram of Figure ^ we plot the mass fraction in 
each phase in Figure H, as a function of the expansion factor 
(normalised to the present value). All phase boundaries are 
fixed, with the exception of the uncollapsed phase - the tem- 
perature threshold is modified to T < {l + z)^K , accounting 
for the adiabatic expansion of the uncollapsed gas. Initially, 
all the gas is in this phase, with a temperature of lOOK . As 
the simulation progresses, some of the gas transfers to the 
shocked phase, as the first virialized haloes become resolved 
{z ~ 10). Gas does not start to enter the galaxy phase (form- 
ing the first objects) until later {z ~ 3). Both the shocked 
and galaxy phases continue to accrete more and more mass 
as the uncollapsed phase becomes progressively depleted. By 
the present day, ~ 18% of the gas is still in the uncollapsed 
phase whilst the shocked gas phase contains about 64% of 
the gas mass and the remaining 18% has cooled into the 
galaxy phase. 

Another important question we have addressed is to 
determine the precise trajectories in the p — T plane, of the 
particles that finally end up in the galaxy phase. Figure ^ 



Figure 4. Two examples of particle trajectories in the p — T 
plane, from the initial epoch (z = 24) to the present day. The 
legends depict their maximum temperatures. The solid path il- 
lustrates the particle that has the highest temperature on ini- 
tial heating (i.e. the first time it crosses the 10*K boundary), 
whereas the dotted path illustrates the evolution of a particle 
which is only heated to lO^K . Marked points indicate redshifts 
z = (10,5,3,2,1,0.5,0.1). 



illustrates the two main paths that the galaxy particles take, 
for the whole duration of the fiducial simulation. The marked 
points indicate the redshifts z — (10,5,3,2,1,0.5,0.1) for 
each trajectory. 

The first example (solid line) shows the expected evolu- 
tion of a galaxy particle - it becomes initially shock heated 
to over lO^K before cooling back down on to the lO^K 
boundary and, on average, progressively increasing in den- 
sity until the final time. The second example (dotted line) 
takes a different path - again, the particle is shock heated, 
but only as far as lO'^K . Both particles undergo a series of 
further heating events and, notably, the second one reaches 
its largest temperature at some later stage when a major 
galaxy merger takes place. 

Quantitatively, we find that only ~ 11% of the parti- 
cles which end up in the galaxy phase were initially heated 
above lO^K for a minimum of one timestep (~ 32% initially 
reached temperatures above 12500K). Most of the galaxy 
mass reaches high density in the smallest haloes close to the 
resolution limit of the simulation. Consequently, their very 
short cooling times causes them to be heated above the lO'^K 
locus and subsequently cool back down within one timestep. 



3.4 Selection Criteria 

To select the "galaxies" , we initially extracted the subset of 
the gas particles that lie in the defined galaxy phase (Fig- 
ure |2|). This was then run through a friends-of-friends group 
finder (Davis et al. 1985), which links particles together to 
form groups until no unlinked particles remain that are 
closer to a group than the linking length, I = bn~^^^, where 
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Figure 5. Projected positions of the fiducial set of 53 galaxies. 
Ttie dots stiow the particles that are linked together by the group 
finder and the circles illustrate the size of the galaxies, centred 
on the median position of each galaxy with radii proportional to 
the cube root of the galaxy mass. 
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Figure 6. Projected positions of the dark matter particles in the 
fiducial simulation selected by the halo finder, Only haloes with 
more than 100 particles are shown. The circles are centred on each 
halo's centre of mass, and their radii are set equal to the distance 
from the centre to the most outlying particle of the halo. There 
are 31 haloes in total in the fiducial catalogue. 



n is the global mean number density. We used the value 
h = 0.1 {I - 30/i"^kpc), which gives a linking length of the 
order of the final S2 softening length. Since the clumps are 
very tight and are well separated, the output of the group 
finder is insensitive to variations in the choice of b: changes of 
the order of a few percent in mass were noted when varying 
6 by a factor of two in each direction. 

Once this catalogue was generated it was truncated by 
throwing away objects containing fewer than A'sph particles. 
Objects below this mass are poorly defined because the SPH 
algorithm does not sample their properties well. For the fidu- 
cial simulation, 16 groups out of 69 were discarded in this 
manner, making up ^ 4% of the mass in the galaxy phase. 
This constitutes most of the residual material in the galaxy 
phase that is not part of the final galaxy catalogue. 

We have also checked the effect of varying our defined 
limits of the galaxy phase. Since we make use of the total 
mass of gas in this phase, it is important to ensure that 
we have selected a region within which most of the parti- 
cles are part of objects in the galaxy catalogue. The most 
important boundary is the lower limit to the gas density 
(p/ (p) = 10'^). Since the 10*K feature is visible down to the 
lowest densities, there are particles that have gone through 
the cooling process that lie outside our galaxy phase. Lower- 
ing the density limit by an order of magnitude caused only 
a 3% increase in the total number of particles grouped into 
galaxies. On the other hand, the fraction of the galaxy phase 
present in the final galaxy catalogue dropped from 96% to 
80%. Reducing the density by the same amount again caused 
no further change in the former quantity but the complete- 
ness again decreased to 56%, indicating that at these low 
densities there is contamination from the diffuse halo gas. 
The temperature boundary for the galaxy phase includes all 
relevant particles for p > 10^ (p), since the hotter particles 



are still cooling down. As a check however, we lowered the 
upper limit on the temperature from lO^K to 12000K . This 
produced a change of less than 1% to the completeness and 
lowered the number of galaxies in the catalogue from 53 to 
52. 

The positions of the fiducial galaxies are illustrated in 
Figure^. Here we plot the projected positions of the particles 
linked together by the group finder. The circles illustrate 
the size of the objects, centred on the median position of 
the linked particles in each group, with radii proportional 
to the cube root of the number of particles linked. 

The dark matter haloe s were identified by u sing a spher- 
ical overdensity algorithm (Lacey & Cole 1994) on the posi- 
tions of the dark matter particles in the fiducial simulation 
endstate. The algorithm makes an estimate of the local den- 
sity which is then used to find the centre of overdense re- 
gions. Successive spheres are placed around this centre with 
growing radii, until the mean internal density equals some 
value - we assume the usual value of 178 times the mean 
density, as predicted by the spherical collapse model for the 
virial radii of haloes in an S7 = 1 universe. Iterations are 
performed, with the centre redefined as the centre of mass 
of the object, until the difference becomes small. Substruc- 
ture within haloes was ignored for the purposes of this paper. 
The halo finder produces a catalogue of all haloes with more 
than Nh = 100 particles in the fiducial case, corresponding 
to a mass of ~ 8 x 10^^h~^ Mq. Since the dark matter dom- 
inates the dynamical evolution, a lower limit in halo mass 
(rather than particle number) was applied when compar- 
ing different simulations. Hence, for simulations with poorer 
mass resolution, the lower limit on the number of particles 
was reduced in order to keep the halo mass threshold con- 
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Figure 7. Phase diagram of all the gas particles within the se- 
lected halo virial radii in the fiducial simulation. Overlaid are the 
boundaries which we have defined in order to split the gas into 
separate phases. 



stant. This limit never went below 1 particles, roughly th e 
lowest number that can be trusted ( Efstathiou et al. 198^ ). 

The output of the halo finder for the fiducial simulation 
is illustrated in figure ^. The circles are centred on the halo 
centre of mass and have radii equal to the virial radius of 
each halo. The dots represent the positions of all the dark 
matter particles within the virial radius of each halo. The 
largest halo has a total mass of ~ 5 x 10^^ Mq, which 
corresponds to a virial temperature of ~ lO^K . 

The p — T distribution of the baryons located within 
the halo virial radii for the fiducial simulation are shown in 
Figure ^ The range of temperatures for the majority of the 
hot phase gas agrees with the corresponding range of virial 
temperatures of the haloes for the mass range of the simu- 
lation (from ~ 5 X 10^ — lO^K). Note that only a negligible 
amount of cold gas falls within the halo boundaries. 



4 SIMULATION COMPARISONS 
4.1 Global Features 

The complete set of simulations that we analysed are listed 
in Table ^, split into three sections. The first section includes 
simulations in which features other than simple parameter 
values were varied. The second section contains simulations 
in which numerical parameters (i.e. those that affect the 
algorithm) were varied. Finally, the third section contains 
simulations in which the physical parameters were altered. 
For most of the parameters, there is a set of three simula- 
tions: the fiducial case and two additional simulations with 
the appropriate parameter varied about the fiducial value. 
The columns indicate the comparison number (13 in total); 
the alteration made to the fiducial parameter set; the frac- 
tion of particles identified in the uncoUapsed, shocked and 
galaxy phases respectively; the number of galaxies in the fi- 



nal catalogue; the number of haloes in the final catalogue; 
and the completeness, which we define as the fraction of par- 
ticles in the galaxy phase that appear in the final object set. 
For discussion, we label the sections general, numerical and 
physical respectively. 

4.1.1 General Comparisons 

The first comparison was designed to probe effects that 
might be present due to the fact that runs were performed 
on different computer architectures. The fiducial simulation 
was run on a Sun Sparcstation, whereas most of the compar- 
ison runs were performed on a Dec Alpha. Hence, we have 
analysed an identical simulation to the fiducial case that was 
executed on the latter platform. Differences arise from the 
system-specific changes made to the code by the compiler. 
Clearly, this change can introduce variations in the phase 
fractions at the level of ~ 1 — 2%. Furthermore, the number 
of galaxies in this comparison differs by 2 and the number 
of haloes differs by 1 . 

Comparison 2 was performed in order to investigate 
the numerical effects introduced by perturbing the particles 
slightly from their initial positions. We set up two simula- 
tions to test this. The first also tested the periodicity of the 
algorithm, by translating the particles by a fixed amount (in 
this case 1 FFT cell) before starting the simulation, then 
moving them back at the final time. The second was per- 
formed in order to examine how small inaccuracies in the 
initial positions affect the results, by adding a random per- 
turbation (tiny compared to the initial displacements) to 
each particle position in the initial conditions. Variations in 
the phase fractions and number of objects are similar to the 
first comparison, hence the two comparisons set a baseline 
for the size of variations that can be introduced solely by 
noise. Simulations that only differ by this amount should be 
viewed as being indistinguishable from the fiducial case. 

For comparison 3, we have added a run that uses a cool- 
ing function in tabulated form, as given by Sutherland & Do- 
pita (1993), rather than a series of power-law fits. The two 
cooling functions used are shown in Figure |l| for the fiducial 
metallicity, Z/Zq = 0.5. The power-law fits the tabulated 
function relatively well for the range of relevant tempera- 
tures (10* — lO^K), although the latter shows a significant 
enhancement in A(r) around lO^K , (as much as a factor of 
3) due to helium line cooling. However, as Table ^ shows, 
there only appears to be a slight enhancement in the amount 
of gas in the galaxy phase. 

The final comparison that was made in this section was 
between the fiducial simul ation and the version of hydra 
used by Pearce et al. ( 1999| ) . The latter contains a change in 
the SPH algorithm t o reduce cooling in large haloes (see also 
Thacker et al. 1998). It is based on a multiphase approach, 



in which the cold, dense particles are decoupled from the hot 
medium: gas above 10"'K does not see gas below 12000K for 
the purposes of calculating the gas density and forces; for 
all other gas particles, they are calculated in the standard 
way. This procedure ensures that infalling cold diffuse par- 
ticles still feel the halo's outer accretion shock and allows 
hot haloes to cool at a rate determined by only the hot gas 
and galaxies to dissipate energy, merge and feel drag as they 
move around the halo environment. The effects of this alter- 
ation on the phase distribution are evident in the table: the 
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Table 2. Details of the set of simulations analysed in this paper, following the various definitions and selection procedures discussed in 
Section ^ The first column is the comparison number given to the simulation, with the second column indicating the alteration made to 
the fiducial initial conditions. The next three columns give the fraction of baryons in each of the three phases defined in the text. Columns 
six and seven give the number of galaxies and haloes respectively that remain in the catalogue after the various selection procedures 
have been applied. The final column is the completeness, i.e. the fraction of baryons in the galaxy phase that are present in the final 
galaxy catalogues. Bold type has been used to denote values that deviate significantly from the fiducial case (two per cent for the phase 
fractions, five percent for the galaxies, haloes and completeness). 



No. 


Alteration 


Uncollapsed Phase 


bhocked Phase 


Galaxy Phase 


Galaxies 


Haloes 


Completeness 





Fiducial Simulation 


0.18 


0.64 


0.18 


53 


31 


0.96 


1 


Architecture 


0.16 


0.65 


0.19 


51 


32 


0.96 


2 


Jr^osition wrtset [^L ceiij 


0.17 


0.65 


0.18 


50 


31 


0.96 


2 


Position Offset (random) 


0.18 


0.64 


0.18 


55 


30 


0.97 


3 


New Cooling Function 


0.18 


0.62 


0.20 


55 


31 


0.98 


4 


Cold gas decoupled 


0.18 


0.67 


0.15 


51 


31 


0.95 


5 


5max — 0.6 


0.17 


0.64 


0.19 


53 


30 


0.96 


6 


K = 0.5 


0.18 


0.63 


0.19 


52 


31 


0.97 


6 


K = 2.0 


0.06 


0.78 


0.18 


53 


31 


0.95 


7 


c I. — ll 

eo = on kpc 


0.16 


0.69 


0.15 


45 


31 


0.95 


7 


eo = 20/i"^kpc 


0.18 


0.63 


0.19 


53 


30 


0.97 


8 


NspH = 16 


0.18 


0.56 


0.26 


97 


31 


0.98 


8 


NsPH = 64 


0.17 


0.70 


0.13 


28 


31 


0.89 


9 


Af = 2 X 16^ 


0.20 


0.80 








29 





9 


A = 2 X 243 


0.17 


0.74 


0.09 


14 


30 


0.71 


9 


A = 2 X 48^ 


0.17 


0.57 


0.26 


141 


34 


0.98 


10 


1 + 2i = 10 


0.17 


0.65 


0.18 


55 


32 


0.97 


10 


1 + 2i = 50 


0.19 


0.63 


0.18 


51 


29 


0.96 


11 


= 0.03 


0.18 


0.69 


0.13 


49 


31 


0.94 


11 


= 0.12 


0.18 


0.58 


0.24 


53 


31 


0.97 


12 


Z = O.OZq 


0.18 


0.78 


0.04 


9 


30 


0.44 


12 


Z = I.OZq 


0.18 


0.60 


0.22 


54 


31 


0.97 


13 


T; = OK 


0.18 


0.64 


0.18 


51 


31 


0.95 


13 


Ti = lO'^K 


0.22 


0.61 


0.17 


50 


30 


0.96 



fraction of gas particles in the galaxy phase is reduced by 
3% (with ~ 2% due to the largest object) while the fraction 
in the shocked phase increases to compensate. 



4-1-2 Varying Numerical Parameters 

The first numerical parameter we varied (comparison 5) 
was the initial (maximum) softening, Smax. Increasing Smax 
moves the changeover from a comoving to a physical soften- 
ing to higher redshift. We doubled the value of the maximum 
softening to Smax ~ 0.6, so that for the same final soften- 
ing length, the changeover now occurred at 2: ~ 3 rather 
than at z ~ 1. More importantly, the softening now has a 
larger comoving value (equivalent to a Plummer softening of 
~ 40 /i~^kpc) before the changeover. The galaxies that form 
before this epoch are more loosely bound (increasing their 
effective volume), and therefore are susceptible to greater 
variations in their mass, due to ram-pressure stripping, tidal 
disruptions or accretion as they move through the hot halo 
gas. No significant changes in the final phase fractions are 
present however. 

The timestep normalisation value, n (comparison 6) was 
varied by a factor of 2 about the fiducial value of 1. This ef- 
fectively doubles or halves the length of each timestep, caus- 
ing the simulation to take roughly twice or half as many 
timesteps to run to z = 0. Doubling the fiducial value in- 
troduces inaccuracies in both the positions and velocities of 
the particles. The phase fractions show that there is signifi- 
cantly more shocked gas in the k = 2 simulation than in the 
fiducial run, at the expense of the uncollapsed gas. The in- 



creased efficiency at which the gas is heated can be explained 
by the fact that the larger timesteps in the k = 2 run in- 
creases the error in particle positions and velocities, which 
leads to larger values in the velocity divergence, increasing 
the efficiency of gas heating from viscous interactions. 

As was hoped, the value k = 0.5 had no significant 
effect on the endstate properties of the galaxies and haloes, 
vindicating CTP95's choice of «: = 1. 

Comparison 7 consisted of halving and doubling the 
value of the final softening length about the fiducial value, 
eo = 10 /i~^kpc. Changing the value of eo only has an effect 
once the physical softening length has dropped below the 
maximum defined by Smax. For eo = 5 h~^kpc this occurs at 
z = 3 whilst for eo = 20 h~^kpc the softening is always co- 
moving. Doubling the fiducial value of eo has no significant 
effect on the final number of galaxies and the fraction of gas 
in each phase. However halving the final softening length 
systematically decreases the number of galaxies by ~ 15%, 
with the residual material still in the shocked phase. This 
might be due to the eff ects of two-body relaxatio n artificially 
heating gas particles (Bteinmetz & White 1997). For a dark 
matter halo with 50 particles that forms ai z — 0, the two- 
body relaxation times are approximately 10.4, 7.5 & 5.9 Gyr 
for eo — 20, 10 &5 /i~^kpc respectively. For eo = 5 /i~^kpc, 
the relaxation time is less than half the age of the universe. 

The value of A^sph was varied in comparison 8. Since 
this parameter sets the resolution at which the density field 
is evaluated, large overdensities are progressively smoothed 
out as the value of A^sph increases. It then becomes increas- 
ingly difficult to form galaxies, since the emissivity is propor- 
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tional to the square of the gas density. On the other hand, 
using a smaller value of A^sph increases the noise arising 
from the discrete nature of the mass distribution. Thus a 
compromise in the value is sought. It is common to assume 
A^SPH = 32, which is our fiducial value, and this comparison 
looks at the effects introduced by changing this number by 
a factor of two. The number of galaxies in each simulation 
clearly shows the expected trend, with more galaxy phase 
gas and more galaxies if A'^sph = 16 than if A'^gpH ~ 64. This 
material comes from the shocked phase as the uncoUapsed 
gas fraction is more or less constant in all cases. The sig- 
nificant drop in completeness for the A^sph = 64 simulation 
is due to 11 groups being discarded because they are be- 
low the imposed resolution threshold (i.e. for these galaxies, 
Ng < Nsph). 

Comparison 9 consisted of changing the mass resolution 
by varying the total number of particles. The initial power 
spectrum was truncated appropriately to account for the 
varying Nyquist frequency of the particle mesh. The four 
runs have N = 2 x (16,24,32,48)^ particles respectively, 
where the 2 indicates the number of species. The masses of 
the two particle species in each of these runs are given in 
Table |^. The three runs with the largest A'^ were selected to 
quantify the result of increasing the resolution of the simula- 
tion. We also added an A'' = 2 x 16^ run to probe the regime 
where artificial two-body heating of the gas particles by the 
dark matter particles should significantly affect the cooling 
rate of the gas. Figure ^ shows the critical dark matter par- 
ticle mass (defined as the mass where the cooling time of 
the gas equals the two-body heating time) as a function of 
gas temperature, for the metallicities and bar yon fractions 
studied in this paper ( ^teinmetz fc White 1997 ). We have as- 
sumed the gas fraction to be equal to the appropriate value 
of Q.h, although the values can be scaled, if desired. The 
Coulomb logarithm is assumed to be 5, appropriate for the 
largest haloes in the simulations studied in this paper (and 
therefore demonstrating the worst case scenario for this ef- 
fect). Overlaid as dotted horizontal lines are the mass of the 
dark matter particles for the simulations in this comparison. 
For Q,b = 0.06 and Z — 0.5Zq, both the fiducial simulation 
and the N = 2x 48'^ simulation are acceptable for the range 
of temperatures relevant to this work, and the N = 2 x 24^ 
simulation is borderline. Since the N — 2 x 16^ run lies well 
above the critical line, the amount of gas cooled in this simu- 
lation should be significantly affected by two-body heating. 
The effect is indeed severe: only 24 particles make it into the 
galaxy phase and no galaxies form at all. (Note that since 
the cooling time at the resolution limit of this simulation is 
^ 10* years - a small fraction of the age of the universe - in 
the absence of artificial heating effects, we should have been 
able to form a population of cooled objects.) 

For the remaining three simulations, the amount of un- 
coUapsed gas remains roughly constant, and the biggest 
change is in the increasing amount of gas that cools from 
the shocked phase to the galaxy phase with larger A'^. This 
is refiected both in the galaxy fractions and in the number 
of galaxies present in the final catalogues. This effect is ex- 
pected when varying the mass resolution of the simulations. 
When A'^ is increased, the sampling of the density field im- 
proves and hence smaller mass haloes are able to be resolved. 
The smallest haloes in the simulation will, on average, form 
first, and therefore have higher density contrasts than the 



Table 3. The masses of each particle species for the runs per- 
formed with varying values of N. 



N 

2 X 16^ 
2 X 24^ 
2 X 32^ 
2 X 48^ 



A/dark(lOlO/l-lAf0) MgaB(lOlO/l-lA/0) 

6.38 0.41 



1.89 
0.80 
0.24 



0.12 
0.05 
0.015 



haloes that form later. As discussed in Section i.3, most of 



the galaxy mass is accumulated in these objects. 

Finally, we investigate the effect of changing the ini- 
tial redshift of the simulation by considering runs with 
1 + Zi = 10,25 and 50 respectively (comparison 10). The 
main potential pitfall when choosing Zi is that choosing too 
small a value will suppress the amount of small-scale power, 
causing objects on small scales to form too late. Hence the 
initial redshift should be high enough that scales of the or- 
der of the initial softening should not have already become 
non-linear (i.e. 5 <g 1). For the 1 + Zi = 50 run, the fidu- 
cial boundaries of the uncoUapsed phase are sufficient to 
segregate the uncoUapsed gas from the rest, while for the 
1 + Zi = 10 run, the temperature boundary is rescaled by a 
factor of 6.25 , i.e. (25/10)^. Evidently, there are no signif- 
icant changes in any of the gas phases or in the number of 
galaxies, when comparing all three runs. 

4-1.3 Physical Comparisons 

The mass fraction of baryons was varied in comparison 11, 
forming the set of values S7b ~ (0.03,0.06,0.12). Varying 
the baryon fraction alters the number density of ions and 
free electrons per simulation gas particle (and therefore the 
masses of both the gas and dark matter particle) , and conse- 
quently changes the cooling rate, since the emissivity scales 
as n^. Again (see Figure ^), it is important that the baryon 
fraction be consistent with the critical mass of the dark mat- 
ter particles, given both the metallicity and the value of A'. 
Clearly for our choice of parameters, the value f2b ~ 0.03 
(Z = 0.5, N = 2x 32^) is borderline and this wiU contribute 
to the reduction in the number of galaxies formed. Increas- 
ing the baryon fraction has little effect upon the number of 
galaxies but they are systematically heavier. 

Comparison 12 was performed in order to examine the 
effect of changing the metallicity of the gas. Each particle 
is assigned a constant amount of metals (i.e. of nuclei that 
have an atomic number greater than 2). The metallicity af- 
fects the shape of the cooling function, as illustrated in fig- 
ure |l|. Also affected is the critical mass of the dark matter 
particles (which depends on the cooling function) as shown 
in Figure ^. The simulation with Z = is clearly on the 
wrong side of the line and the galaxy phase fractions are 
consequently affected: slightly more material cools with the 
higher metallicity, but the run with no metals shows a catas- 
trophic drop in the amount of material in the galaxy phase, 
the number of galaxies and the completeness. 

Finally, in comparison 13, we created a set of runs with 
different initial gas temperatures in order to measure how 
the dynamical and thermal evolution of the gas depends on 
its initial thermal state. We included two runs, one with 
Ti = OK and the other with Ti = lO'^K . AUhough these pa- 
rameter choices may seem hard to justify in physical terms. 



© 0000 RAS, MNRAS 000, 000-000 



Parameter Tests Within Cosmological Simulations of Galaxy Formation 11 



-16 






















L n^ = 0.03 _ 


i ; 

nj^:=o.o6 _ 


/ 



1 1 



10 



o 



6 

log (T/K) 



6 

log (T/K) 



6 

log (T/K) 



Figure 8. The critical dark matter particle mass for spurious 2-body heating to exceed cooling, as a function of temperature ( Steinmetz 
& Whit! 1997). The three solid curves in each plot correspond to different gas metallicities (from top to bottom, Z/Zq = 1.0,0.5,0.0, 



with the fiducial metallicity, 0.5, in bold). The different panels assume a local gas fraction, /, corresponding to the values of the global 
baryon fraction analysed in this paper, = 0.03, 0.06, 0.12. The dotted horizontal lines give the masses of the dark matter particles for 
N = 2x X'-^, where X is the label on each line. A value of 5 was assumed for the Coulomb logarithm. 



we decided to stick with our approach of varying a single pa- 
rameter at a time. The fiducial boundaries of the uncollapsed 
phase are sufficient for the run with Ti = OK , however the 
boundaries are redrawn for the run with Ti = lO'^K , such 
that s < 2.5 and p/ (p) < 50, where s = logi(,(r/(p/ {p)f'^) 
measures the specific entropy of the gas. (Using these con- 
straints makes it easier to pick out the uncollapsed gas since 
it is isentropic.) For the run with Ti = OK , there are no 
significant differences in the measured quantities, however 
the run with T — lO^K shows a small increase in the uncol- 
lapsed gas phase at the expense of the shocked gas. For this 
run, the gas cools as expected from adiabatic expansion until 
2 ~ 5. However by z = 3 the uncollapsed phase is at tem- 
peratures below lO^K , over an order of magnitude below 
the adiabatic temperature. At these redshifts, the density 
and the cooling rate of the gas is high enough to allow it 
to radiate a significant fraction of its energy away. Conse- 
quently, only a small fraction of the gas is too hot to undergo 
gravitational collapse. 

4.2 The Mass Distribution of Galaxies 

In Figure ^ we plot the cumulative mass functions of galaxies 
(i.e. the number density of galaxies greater than a specified 
baryon mass, Mg), with the fiducial simulation always plot- 
ted as the solid line. The dynamic range for these simulations 
is about 2 orders of magnitude in abundance and about 1.5 
orders of magnitude in mass. The functions all have a char- 
acteristic shape: a shallow slope at low masses, turning over 



to a steeper slope at higher masses and a large tail at the 
high mass end. The latter feature is simply due to the exis- 
tence of a single massive galaxy that has a baryonic mass of 
~ 5 X lQ^^h~^ Mq. Galaxies this massive are extremely rare 
in the universe; we shall discuss this object below. 

The biggest differences in the mass functions occur in 
the physical comparisons, in which the gas cooling rate 
varies. The simulation with Z = Q fails to produce as many 
galaxies as the fiducial simulation - the abundance is down 
by a factor of 5 and the galaxies that do form are much less 
massive. The large differences seen in the runs with varying 
baryon fraction are mainly due to the mass of the individ- 
ual gas particles being different by a factor of two. However, 
when the galaxy masses are scaled to take out this effect, 
there is still a systematic difference between each function. 
This is due to the effect the baryon fraction has on the cool- 
ing rate of the gas. 

Varying the number of SPH neighbours shows a pro- 
gression in the smallest resolved galaxy mass due to the 
selection criterion that the number of constituent particles 
Ng > NspH- There are fewer galaxies for larger A''sph and 
they are systematically less massive. This is a refiection of 
the fact that the algorithm has to smooth over a larger range, 
which consequently leads to a poorer resolution of the den- 
sity field. Since emissivity scales as n^, the lower cooling rate 
inhibits the galaxy formation process - an entirely numeri- 
cal effect. Varying the value of N shows similar results for 
the same reason, i.e. a larger number of particles enables the 
SPH algorithm to resolve higher gas densities. 
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Figure 9. Cumulative galaxy mass functions for the full set of simulation comparisons. The plot gives the comoving number density of 
galaxies with baryon mass greater or equal to Alg. The fiducial simulation is always plotted with a solid line-type. 



The remaining differences are relatively small. Chang- 
ing the cooling function to the tabulated form of Sutherland 
& Dopita (1993) produces more cooling around 10^ K , boost- 



ing the cooling rate and systematically increasing the mass 



of the galaxies by ~ 10%. However, there is no significant 
difference in the shape of the mass function. 

Comparison 4 (in which the cold gas is decoupled from 
the hot gas) also shows a small difference in the mass func- 
tions. This difference is dominated by the largest object in 
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the simulation volume, which is reduced in mass by a factor 
of 3 when the cold, dense gas is decoupled from the hot gas. 
Notably, this alteration to the algorithm has a negligible 
effect on the smaller objects. 

The comparison set in which the final softening length 
is varied shows an offset for eo = 5/i~^kpc. This simulation 
produces galaxies that are systematically lighter than the 
simulations with larger values of eo, as expected if they are 
being affected by two-body heating from the dark matter 
particles. 

4.3 Comparison of Matched Galaxies and Haloes 

To compare the individual properties of the galaxies and 
haloes, we performed object-to-object matching between 
each simulation and the fiducial run. Since the number of 
objects in each catalogue is different in nearly every case, 
there is not a straightforward one-to-one correspondence 
between the sets. To circumvent this, we always took the 
smallest object set first and compared the larger one to it. 
When more than one object in the larger set was matched 
with the same object in the smaller set, the closest pair was 
selected. Generally, this was not a problem since the dis- 
carded objects are usually at considerably larger distances 
than their matched object. However, this method will in- 
evitably produce outliers in both separation and mass dif- 
ferences due to the fact that, on occasion, an object in one 
simulation is fragmented into smaller objects in another. 
This prompted us to use the median and semi-interquartile 
ranges as statistical measures of the offset and scatter of 
any relative quantity we measured. Positions of the galaxies 
were defined as the median position of the linked particles 
and, for the haloes, their centre of mass. Below we detail the 
individual comparisons we made. 

^.3.1 Galaxy-galaxy displacements 

The first property we examined was the scatter in the po- 
sitions of the matched galaxies. The results are illustrated 
in Figures and |ll|, which show the displacement between 
each matched pair (in /i~^kpc) plotted against the average 
object mass. The solid line illustrates the median separation 
for each comparison. The number of galaxies matched and 
the median plus semi-interquartile ranges of the displace- 
ments are presented in Table ^, columns 3, 4 & 5. It is clear 
from the comparison between the simulations run on differ- 
ent computer architectures that there is an intrinsic offset 
in matched positions at the level of a few softening lengths. 
Furthermore, a significant number of matches exceed this 
level by as much as a factor of four. As a consequence of 
this, we regard as significant only those comparisons that 
have median offsets in excess of the level seen for these two 
simulations. 

The first set to show a significant change is the simula- 
tions with different values of A''. Varying the mass resolution 
causes differences in the number of objects formed and also 
affects the background dark matter distribution, so it is of no 
surprise that the displacements show a significantly larger 
amount of scatter in these comparisons. Varying the value of 
A^SPH also produces very different numbers of galaxies but 
this does not affect the median separation at all, demon- 



strating that it is the difference in the dark matter that pri- 
marily drives the increased scatter in the galaxy positions 
rather than mergers or interactions between the different 
numbers of galaxies. 

A systematic difference in the matched positions is also 
introduced when the initial redshift of the simulation is 
changed: the median displacement increases by a factor of 
two in both comparisons. Variations in the initial positions 
and velocities of the particles leads to asynchrony in the 
subsequent spatial trajectories of the galaxies. 

Regarding the physical comparisons, substantially in- 
creasing the initial temperature of the gas leads to a sig- 
nificant effect. Particles that have been supplied with this 
much thermal energy naturally have greater pressure sup- 
port, which directly affects the equation of motion, and 
therefore the particle trajectories. 

4-3.2 Galaxy-galaxy masses 

We have also analysed the scatter in the masses of the 
matched galaxy pairs. Figures ^ and ^ illustrate the re- 
sults by showing the logarithm of the mass ratio of each 
matched pair plotted against their average baryonic mass, 
always defined as the variant relative to the fiducial case. 
The solid line illustrates a ratio of unity (i.e. equal masses) 
and the dashed line illustrates the median value of the mass 
ratio. The median and semi-interquartile ranges are given 
in columns 6 & 7 of Table ^ 

Adopting the tabulated cooling function causes an in- 
crease in the mass of each galaxy, because of an overall en- 
hancement in the cooling rate around the temperatures most 
appropriate for these simulations (~ 10^ — lO^K). Decou- 
pling the hot gas from the cold gas causes the masses of 
the galaxies to decrease slightly with only a few exceptions, 
notably the largest object. Decreasing the softening length 
from eo = 10 to 5/i~^kpc also has the effect of decreasing 
the masses of the galaxies. Similarly, a larger value of A^'sph 
causes the masses of the galaxies to be systematically lower, 
due to its effect on the density field and hence the cooling 
rate. 

The dispersion in the mass ratio for the comparisons 
with varying A*' is large because of the large scatter in object 
positions for these runs. This makes spurious matches much 
more likely. For A = 2 x 24^ the objects are significantly 
less massive, whilst for A = 2 x 48^ the objects are only a 
little more massive than the fiducial simulation. 

The largest differences in galaxy mass result from al- 
tering the cooling rate via the physical parameters, fib and 
Z. For the comparisons with varying baryon fraction, the 
galaxies are more than a factor of two heavier for the larger 
value of fib- Furthermore, there appears to be a slight trend 
for increasing mass excess for heavier galaxies but this is dif- 
ficult to measure due to the small dynamic range in mass. 
Varying the metallicity produces a similar effect to varying 
fib. 

4-3.3 Dark matter halo-halo masses 

The method used for comparing masses of matched galaxy 
pairs was also performed on the dark matter haloes. 
Columns 8 & 9 in Table quantify the differences. Over- 
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Table 4. Statistical measures of the scatter present in all simulation comparisons of the matched galaxies and haloes. The first two 
columns detail the particular comparison and the third column lists the number of galaxies in each catalogue that were matched with 
the fiducial simulation, using the method detailed in the text. The rest of the columns list the median and semi-interquartile range (siqr) 
for the measured scatter in galaxy separation (in h~^kpc), galaxy mass, halo dark matter mass and halo baryon mass respectively. 
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all, very tight correlations are seen for most of the simu- 
lations, as the dark matter is the dominant source of the 
gravitational potential and so is not significantly affected by 
changes in the gas dynamics. The simulation with fib = 0.12 
shows the largest difference - here the contribution from the 
gas starts to make a difference to the properties of the dark 
matter haloes. 

4.4 The Distribution of Baryons in Haloes 

The final item we investigate in this study is how parameter 
variations affect the final distribution of the baryons within 
haloes, as a function of the halo virial mass. Specifically, we 
looked at the mass fraction of baryons present in the galaxy 
phase, the ratio of baryon masses for matched haloes and the 
total halo baryon fraction normalised to the global value. 

4-4- i Fraction of baryons in galaxies 

Figure |l^ illustrates the fraction of baryons in the galaxy 
phase as a function of the halo virial mass, for each com- 
parison. Because of our small samples, we have binned the 
data as histograms; the last bin corresponds to the single 
largest object. This statistic is remarkably stable; the first 
two comparisons show virtually no difference. The trend is 
for higher mass haloes to have less of their mass in galax- 
ies. The mass fraction varies from ~ 0.65 in the smallest 
haloes down to ~ 0.2 in the highest mass halo. The compar- 
ison between the fiducial simulation and the case in which 
the galaxy phase is decoupled from the hot gas reflects the 
reduction in the total amount of galaxy phase baryons in 
the latter - about a factor of three for the galaxy in the 
largest halo (~ 5 x IO^^/i'^Mq), but only ~ 10% for haloes 
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with masses of ~ 10^^ Mq. Other comparisons that give 
significant changes in the mass fractions show nothing un- 
expected. 



4.4-S Local halo baryon fraction 

The variations in the total mass of baryons within the virial 
radius of matched haloes are quantified in columns 10 & 11 
of Table ^. The median differences are slightly higher than 
for the dark matter, but smaller than the galaxies alone. In- 
terestingly, the sign of the median difference is not always 
the same as the equivalent comparison for the galaxies, indi- 
cating that fluctuations in the total baryon content are not 
due to the galaxies alone. Significant differences are most ev- 
ident in the simulations in which parameter changes have af- 
fected the cooling rate. Inducing a higher cooling rate causes 
the gas to lose some of its pressure support and therefore to 
become more compressed within its halo. The simulation 
with eo = 5/i~^kpc shows a significant change in the same 
direction as the change in the galaxy masses, supporting the 
idea that the gas is being more strongly heated from two- 
body encounters with the dark matter particles. 

We also analysed the total baryon fraction in haloes as 
a function of the halo virial mass. The results are displayed 
in Figure |l^. The fraction is normalised to the global value 
of ilb for each simulation. The fiducial simulation has halo 
baryon fractions ranging from ~ 93% of the global value in 
the smallest haloes to ~ 82% in the largest. All simulations 
agree with the fiducial values to within ~ 20%. Notably, 
not all simulations show a monotonically decreasing baryon 
fraction with halo mass. However, the scatter is considerable 
because of the small number of objects in each bin (with 
only one in the highest mass bin). The differences reflect 
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Figure 10. Separations between matched pairs of galaxies in the fiducial simulation and all other simulations, plotted against the average 
galaxy mass. The labels indicate the parameter that varies in each plot. The median scatter is shown as a solid horizontal line. Matches 
above 100/i~^kpc are not plotted. 



the varying amounts of material that has cooled. Note, for 
example, the drop in the first bin for Z = 0.0, for which 
galaxies do not form in small haloes and the rise in the first 
bin for Ai"spii = 16. For almost all of the simulations and 
across the entire range of halo masses, the baryon fraction 
is less than the global value. The baryon fraction in the 
largest halo in this simulation (0.82 for the fiducial case) is 
somewhat smaller than the values, 0.9-1.0, found in the rich 
galaxy cluster simulations presented by Frcnk ot al. (1999). 
The value for our largest cluster, however, is noisy as indi- 



cated by the first comparison in the figure and it depends on 
metallicity, as indicated also in the figure. Furthermore, this 
cluster is approximately 10 times less meissive than the clus- 
ter in Frciik ct al. and, as shown by Poarco et al. (in prepa- 
ration), there is a weak trend of increasing baryon fraction 
with mass for clusters in the SCDM cosmology. 
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Figure 11. Separations between matched pairs of galaxies (continued). 



5 SUMMARY OF THE COMPARISON 
RESULTS 

We now summarise the main effects seen in each comparison, 
in turn. 

In comparison 1, we looked at the effects introduced by 
running the simulation on different computer architectures. 
No significant differences were found between the two simu- 
lations. However, the displacements between matched pairs 
of galaxies at z = were found to be of the order of an S2 
softening length. It is therefore incorrect to trust either the 



positions of individual galaxies or their trajectories on scales 
smaller than this. 



Comparison 2 was intended to look at the effects of 
small displacements in the positions of the particles at the 
start of the simulation. Wc ran one simulation with the par- 
ticles initially displaced by 1 FFT cell to test the periodic 
nature of the algorithm. Wo also ran a simulation with tiny 
random perturbations applied to the initial particle posi- 
tions, to test the level at which differences can grow over 
the course of a few thousand timesteps (the typical length 
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Figure 12. The logarithm of the ratio of the masses of matched pairs of galaxies plotted against the logarithm of the mass of the fiducial 
galaxy. Labels indicate the parameter corresponding to each comparison (with the relative difference defined as the second minus the 
first value). The solid line corresponds to equal masses and the dashed line illustrates the median value for each comparison. 



of each simulation). Again, no significant changes were de- 
tected amongst the three simulations. 

We next changed the cooling function from the series 
of power-law fits, used in the fiducial simula tion, t o a set of 
tabulated values from Sutherland & Dopita (1993). A slight 
enhancement in the galaxy masses, at the 10% level, was 
produced in the latter case, although no significant changes 
were produced in the spatial distribution of galaxies. 

In comparison 4, we compared the fiducial simulation 
to one in which the galaxy phase gas is partially decoupled 



from hot gas in the SPH algorithm (see Pearce et al. 1999). 
This modification provides a better estimate of the hot gas 
density near galaxies, specifically it prevents an overestimate 
which leads to an artificial enhancement in the cooling rate. 
Only the largest galaxy in the simulation was significantly 
affected, with its mass being reduced by a factor of three. 

Comparison 5 consisted of doubling the initial (maxi- 
mum) softening length. This results in a period of time in 
which the softening is larger than in the fiducial simulation. 
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Figure 13. The logarithm of the ratio of the masses of matched pairs of galaxies (continued). 



This led to slight differences in the spatial distribution of the 
galaxies, although it left the mass distribution unaffected. 

We altered the size of the timesteps by a factor of two 
each way, in comparison 6, by changing the constant, k. 
Increasing the timestep resulted in much more hot gas at 
the expense of the cold phase. This illustrates the effect of 
the tirncstcp on the dynamics of the particles: introducing 
significant errors into the positions and velocities leads to 
enhanced heating of the gas. The galaxies themselves re- 
mained largely unaffected. No significant change was evident 



when halving the value of k, thus justifying the value /t = 1, 

adopted as fiducial. 

Comparison 7 consisted of changing the final value 
of the softening length, eo- In the simulation with eo = 
5/i^^kpc, the galaxies were about 20% less massive than 
in the fiducial simulation (eo = 10/i^^kpc), possibly due to 
discreteness effects caused by harder two-body encounters. 
A larger value of eo = 20ft~^kpc, which implies a constant 
comoving softening, produced no significant differences. 

In comparison 8, we varied the value of iVgPH, malt- 
ing it larger and smaller by a factor of two, thus looking at 
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Figure 14. The mass fraction of baryons in the galaxy phase plotted as a function of the total halo mass. The fiducial simulation is the 
solid histogram. 



the effects produced when the SPH algorithm uses a differ- 
ent number of neighbours when calculating the smoothing 
length for each gas particle. These changes had a dramatic 
effect on the simulations. The value of iVspH determines the 
minimum mass for which the gas density can be resolved: 
smaller values will resolve smaller masses at the expense of 



larger random fluctuations. We found that the number of 
galaxies increases strongly with smaller iVspH, as did the 
masses of the matched galaxy pairs, although without sig- 
nificantly affecting their positions. 

Comparison 9 consisted of changing the number of par- 
ticles in the simulation, thus altering the mass resolution. 
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Figure 15. The baryon fraction in haloes as a function of the total halo mass, in units of the global baryon fraction. The fiducial 
simulation is the solid histogram. 



To demonstrate the effect of spurious heating by dark mat- 
ter particles, we ran a simulation with a dark matter parti- 
cle mas s abo ve the critical value calculated by Steinmetz & 
White (1997). In this case, gas is unable to cool and form 
galaxies. For the main set of comparisons, we had a suite 
of runs with N ^ 2 X [24,32,48]^. The number and mass 



of galaxies increased with A'^. Improving the mass resolution 
also affected the final distribution of dark matter. The two 
simulations with the highest resolution were in better agree- 
ment than the two with the lowest resolution, albeit that the 
discrepancies in the galaxy phase material (e.g. the fraction 
of cooled gas) were still large. 
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In comparison 10, the initial redshift of the simulations 
was changed (1 + = 10,25,50). The masses of the ob- 
jects were not significantly altered by this change but small 
changes in the spatial distribution of the galaxies were pro- 
duced. 

Comparison 11 involved changing the value of the 
baryon fraction (S7b=0. 03, 0.06, 0.12). The higher value of fib 
causes the cooling rate to increase (since emissivity is pro- 
portional to r^b)- The mass of each galaxy also increases, 
since m is proportional to Q,b. The run with f2b = 0.03 has 
dark matter masses of the order of the critical mass for 2- 
body heating, although the reduction in the cooling rate 
produced by this effect is confused by the fact that the cool- 
ing rate is significantly less anyway. Overall, large differences 
were noted in all analyses in which the baryon fraction was 
varied. 

In comparison 12, we changed the value of the (fixed) 
gas metallicity {Z — 0.0, 0.5, 1.0 in solar units). The run with 
no metals differed dramatically from the fiducial simulation, 
producing very few objects by 2: = 0. This simulation had 
a dark matter mass well above the critical value for 2-body 
heating. The run with Z — I.OZq showed an enhancement 
in galaxy masses relative to the fiducial simulation, although 
the change was far less severe. 

Finally, in comparison 13 the initial gas temperature 
was changed from lOOK to OK and lO^K respectively. No 
significant changes were measured in the former two. Al- 
though the uncoUapsed gas in the latter simulation is on a 
higher adiabat (T ~ lO'^K), this had very little effect on 
the final state of the gas in galaxies and haloes, with the 
exception of introducing significant displacements between 
matched galaxy pairs. 



6 CONCLUSIONS 

In this paper, we have explored how variations in simulation 
parameters and physical assumptions affect the population 
of galaxies which form in cosmological simulations. We used 
the AP''m algorithm to calculate the gravitational potential 
field and SPH with radiative cooling to model gas dynamics. 
We started with a carefully chosen fiducial simulation and 
calculated a number of variants, each differing in only one 
parameter. Galaxy and halo catalogues were selected in a 
consistent manner before analysing the effects that the pa- 
rameter changes had on the thermodynamic state of the gas, 
on the mass and spatial distribution of the galaxies, and on 
the baryon content of haloes. 

By carefully tracking the phase evolution of gas parti- 
cles we have shown that most galaxy material (around 90% 
in the fiducial simulation) originally cools in small haloes 
where the virial temperature is ~ lO^K . As Figure shows, 
the percentage of halo baryons that are cold rises as the halo 
mass falls. Only a small fraction of the baryons cool at later 
stages of the evolution. Mergers of these small fragments can 
reheat the gas briefly as massive galaxies are built up. This 
implies that the most dynamically interesting part of the 
process occurs at the resolution threshold of the simulation 



for most galaxies, but can be much higher in some cases. 
Individual galaxy trajectories should therefore be treated 
with caution. 

The length of timestep recommended by Couchman, 
Thomas & Pearce (1995) is adequate. Halving it caused no 
significant change in the galaxy catalogues, but choosing a 
larger value led to a significant increase in the shocked gas 
relative to the uncoUapsed gas. This can be attributed to 
errors in the particle positions and velocities, resulting in an 
artificial enhancement in the thermal heating of the gas. 

We confirm that simulations with dark matter particle 
masses greater than the Steinmetz & White (1997) critical 



mass for 2-body heating do indeed show a catastrophic re- 
duction in the final amount of cooled gas and in the number 
of galaxies, even when the shortest cooling time is a small 
fraction of the age of the universe. The choice of dark mat- 
ter particle mass should be comfortably below this critical 
value. The mass in galaxies is also strongly affected by the 
choice of baryon fraction and metallicity. 

Pre-heating of the gas to a high initial temperature 
(lO^K at 2 = 24) produces small changes in the final posi- 
tions of the galaxies but has little effect on their final mass. 
A significant fraction of the thermal energy of the gas is ra- 
diated away around z = 3 — 5, allowing the majority of the 
baryons to undergo gravitational collapse. 

Increasing the mass resolution by using a larger num- 
ber of particles causes both the number of galaxies and their 
masses to increase. Improved resolution of the density field 
enhances the cooling rate within haloes, particularly at high 
redshift when the first objects form. As the mass resolution 
of the simulation increases, modelling feedback processes 
correctly becomes critical in order to avoid overcooling of 
baryonic material. 

Lowering the number of neighbours used by the SPH 
algorithm to calculate hydrodynamical quantities has a sim- 
ilar effect on the amount of cooled material to increasing the 
total number of particles in the simulation. Smaller values of 
A'sPH enable higher overdensities to be resolved (therefore 
enhancing the cooling rate in such regions), at the expense 
of more contamination from sampling noise. 

The fraction of the baryonic material that cools within 
a halo and forms galaxies is a very stable quantity which 
varies inversely with halo mass. Reducing the total amount 
of gas affects all scales in a similar way. The ratio of baryonic 
to dark material within the virial radius of haloes is found 
to be around 0.8 ±0.1 of the global value, for the range of 
halo masses studied here. 

In this study we have used four different prescriptions 
for the gravitational softening. We have varied the initial 
comoving softening, from 20/i^^kpc to 40/i~^kpc, and we 
have also used a range of final physical softenings, from 
5^~^kpc to 10/i~^kpc and 20/i~^kpc. With the exception 
of the 5^~^kpc case, this variation had very little effect on 
the galaxy population. The smaller final softening system- 
atically reduced both the mass of galaxies and the mass of 
baryons within haloes by ~ 20%. This effect is probably due 
to the artificial heating of the gas from two-body encounters 
with the dark matter particles discussed above (Steinmetz 



and is therefore sensitive to many numerical and physical 
parameters. 

Numerical noise leads to a scatter in the final galaxy 
positions that is of the order of the (spline) softening length 



White 1997). The fiducial softening is borderline for this 



resolution. 

Finally, in almost every simulation we obtained an 
overly massive galaxy at the centre of the largest halo. With 
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a baryonic mass of ~ 5 x 10^^ A'Iq , such an object is not 
expected within such a small volume of the universe. De- 
coupling the central dense gas from the hot halo material 
effectively reduces t he mass of this ob ject to a more reason- 
able value (see also Pearce et al. 1999). 

We conclude that varying most of the numerical or 
physical parameters required in simulations of galaxy for- 
mation within reasonable bounds produces only relatively 
small changes in the abundance, mass and spatial distribu- 
tion of the galaxies that form. The most important numer- 
ical parameter is mass resolution which has an important 
effect because most of the gas cools in small mass haloes (in 
which the cooling times are shortest). The most important 
physical parameters are those that control the cooling time 
of the gas, the density of baryons and its metallicity. Vary- 
ing these within plausible ranges can substantially change 
the abundance and mass of galaxies. Careful modelling of 
physical processes occurring in small mass haloes where the 
thermodynamic state of the gas is likely to be affected, for 
example, by energy injected by supernovae and stellar winds 
is the next great challenge facing simulators of galaxy for- 
mation. 
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